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ABSTRACT 


A modal approximation for the self and mutual radiation 
impedances has been derived for arrays of spherical 
transducers with small ka values, where ka is the acoustic 
wave.number multiplied by the radius of the sphere. I term 
this the "Modal Pritchard" approximation, as it is related 
to the so-called Pritchard approximation, often employed to 
calculate mutual radiation iitpedances. I investigated the 
utility of the approximate mutual radiation impedance 
expression for three two-body problems (monopoles, aligned 
dipoles, and aligned-linear quadrupoles). ^ For these cases, 
approximate values were found to be in good agreement with 
those obtained using a full Spherical Addition Theorem 
calculation, and are an improvement over the simple 
Pritchard approximation. Additionally, I investigated the 
mutual radiation impedance expression in one particular 
three-body problem. Because of the Modal Pritchard 
approximation's inability to correctly handle scattering, we 
recommend using the full Spherical Addition Theorem 
calculation when scattering is important. 

Finally, I investigated the use of a new finite element 
mesh to calculate the T-matrix for a given transducer. The 
T-matrix relates the incident and scattered waves for a 
single transducer, in an orthogonal (spherical harmonic) 
basis set. The monopole element showed an increase in 
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error, while we saw some improvements in the higher-order 
diagonal elements. Off-diagonal elements, which should be 
zero for a spherical scatter, were satisfactorily small in 
most cases. Although the results were less than favorable, 
I was able to streamline the T-matrix calculation while 
providing a new method of examining the off-diagonal 
elements. 
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I. INTRODUCTION 

The work described in this thesis is part of an ongoing 
effort to irtprove our ability to accurately model densely 
packed underwater acoustic arrays. When an acoustic array's 
operating characteristics are closely predicted we say the 
array behaves well. This behavior is influenced by element 
interaction, characterized by a parameter called mutual 
radiation irrpedance (Ref. 1). A worst case example of "bad" 
behavior is that a transducer element may actually absorb 
acoustic power, which is then absorbed back into the driving 
amplifier, thus possibly causing an overload and electrical 
failure (Ref. 2) . Other examples of bad behavior include 
changes in the radiated source level and deviations of the 
beampattem, steering direction, and side lobe amplitudes. 

Element-to-elemeht interactions are more significant in 
arrays that are packed into a volume of small size, which is 
the trend in low frequency active arrays. We will analyze 
the mutual radiation impedance for spherical sources of 
small ka values, where k is the wave n;imber and a is the 
radius of the sphere. The mutual radiation impedance is an 
important parameter used to calculate the pressure on one 
array transducer due to the radiation and scattering from 
other transducers in the array. 

Additionally, I will use a finite element code to solve 
for the radiating and scattering properties of a single 
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transducer. Therefore, this thesis consists of two 
portions. 

The first part is devoted to deriving a general 
expression for the mutual radiation impedance between two 
identical radiating hard spheres in an otherwise free field 
environment using a normal mode approach. Additionally, an 
expression for the self-radiation impedance, which is the 
radiation loading on a single element in a free field, is 
found. 

I have compared our modal approximation for the mutual 
radiation impedance with the commonly-used so-called 
“Pritchard" approximation for the particular case ka = 1, 
and will show the improvement possible using our modal 
approximation. This improvement is particularly noticeable 
as the center-to-center distance between sources is reduced. 
In this "Simple Pritchard" approximation, as I refer to it, 
the real part of the self-radiation impedance is multiplied 
by a spherical Hankel function of zeroth order to find the 
mutual radiation impedance. I will show how this "simple" 
approximation starts to break down as the transducers are 
spatially brought closer together. 

The Addition Theorem (Ref. 3) provides one the ability 
to mathematically express the pressure radiated from one 
body referenced to the coordinate system of another body. 
Using the Addition Theorem, it is possible to very 
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accurately compute the mutual radiation impedance of a 
sphere in an arbitrarily densely packed array of any number 
of elements. I have used the Addition Theorem as a 
computational standard. 

I also applied our modal approximation to a specified 
three-body problem to further investigate our modal 
approximation, which I will hereafter refer to as the "Modal 
Pritchard" approximation. In this problem we note the 
inability of our Modal Pritchard approximation to handle the 
effects of scattering. Because of this breakdown we 
recommend using the full Spherical Addition Theorem for this 
problem and other cases when scattering becomes important. 
Methods and computer codes for this portion of the thesis 
are provided in Appendixes A and B. 

The second part of this thesis makes extensive use of a 
finite element code, named ATILA, which was written by 
engineers at the Institut Superieur d'Electronique du Nord 
(ISBN). ATILA has been specifically developed to aid in the 
design of sonar transducers, but is also used to model a 
variety of acoustic problems. 

The final portion of this thesis addresses the T-matrix 
method (Ref. 3). In order to model certain aspects of array 
operation such as near-field pressure, it is necessary to 
solve for the radiating and scattering properties of a 
single transducer. We accomplish this calculation by using 
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the ATILA finite element code. The T-Matrix (or Transition 
Matrix) characterizes the scattering properties of a 
specific body by translating incoming pressure waves to 
outgoing ones. 

I will provide a brief description of ATILA's method of 
calculation that distinguishes it from other finite element 
codes. Readers less familiar with finite element codes may 
want to seek further information from reference books such 
as Kwon's "The Finite Element Method using MATLAB" (Ref. 4). 

Previous work with ATILA to calculate the T-matrix for 
a thin-spherical shell has been shown to be in error when 
compared with a verified analytical solution (Ref. 5). 
Various limitations of the ATILA code may ej^lain the errors 
present. I investigated whether improvements are obtained 
by making some corrections tp ATILA's code and by the 
implementation of a new, refined, finite element mesh. This 
new mesh increased the errors to the monopole element but 
did provide some improvement to other T-matrix elements. 
Further irrprovements may be realized with the advent of more 
sophisticated radiation boundary operators. An early 
delivery of an addition to the ATILA code, EQI, which uses 
using a boundary element method to handle the fluid, did not 
occur in time for use with this thesis. Investigation of 
the improvement obtained using ATILA coupled with EQI will 
have to be the subject of future research. 
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In addition to providing some inprovement in the 
calculation of the T-matrix, I have been able to streamline 
the procedure. Output scattered pressures from ATILA were 
post-processed using easy to read MATLAB program scripts. 
The resulting T-matrix calculations were quickly compared 
with the analytical results (also done in MATLAB) . Listed 
in Appendix C are the MATLAB codes used to calculate the 
theoretical values of the T-matrix elements and the T-matrix 
elements calculated from ATILA's scattered pressure results. 

In summary, by deriving an improved analytical 
formulation for the mutual radiation impedance, one can 
quickly investigate the effects of parameter variations, 
such as the distance between elements in the array. 
Additionally since this formulation involves angular 
dependence, I have devised a means of investigating the 
effects of altering the orientation between elements in an 
array. We recommend using the Modal Pritchard approximation 
for approximate solutions when scattering is unimportant, 
and the full spherical Addition Theorem otherwise. 
Streamlining the T-matrix calculation for a given transducer 
will allow follow-on work to be accomplished with less 
startup time. 
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IZ. THEORETICAL BACKGROUND SECTION 

This chapter addresses the development of a modal 
approximation for the self and mutual radiation impedances 
for two spheres in an otherwise free field environment. 


A. BOUNDARY CONDITIONS 

Consider a space and time dependent velocity potential, 
defined as lP’(r, t) = W{r)e^‘^, where the vector r = (r,0,(p) 

represents the spatial dependence. A relationship between 
the pressure (p) and the velocity potential ( *F) states p = 
-pjd'i'/dt (Ref. 6), so that p = -jOJp^W. 

Taking the gradient and then the dot product with the 

radial unit vector, we obtain dW/dx - -II (japJ dp/dr. By 

} ' ' 

noting that the velocity vector is defined to equal the 
gradient of the velocity potential (u = , we can derive 

a boundary condition for a surface normal velocity. For a 
sphere of radius a, for example, evaluation of dW/dx - - 
1/ {jCDpJ ^f0x at r=a provides an expression for the radial 
velocity on the sphere's surface. This will be a boundary 
condition for our analysis. 



j 

oypo 
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B. PROBLEM SETUP AMD ADDITION THEOREM 

Consider two identical spheres (radii = a) spaced with 

a center-to-center distance = d. Letting and 

be the spherical coordinates relative to the 
centers of spheres one and two, respectively. With the two 
spheres radiating in distinct normal modes, the radial 
velocities are expressed below. 

Where ^ modal eigenfiinction, p“ is the 

Legendre function, and are modal amplitudes. 

The pressure from each sphere can be mathematically 
expressed as shown in Equation (3) (Ref. 3) . Since we are 
assximing a e^'^ time dependence, we use a spherical Hankel 
function of the second kind. If we were to assiume a time 
dependence of then we would instead use a spherical 

Hankel function of the first kind. After Equation (3), we 
drop the Hankel function's (2) superscript. 



Given the Addition Theorem, we can express the pressure 
from sphere two relative to sphere one (see Equation (4), 
Ref. 7). The appendix within Ref. 7 provides the definition 
for the function a{v,p, t,fi, s) . The angles 0^^ and are 
formed by describing the location of sphere one's center 
relative to sphere two's center (see Figure 1). The prime 
in the last summation means the index is incremented by two. 

t i 't 

v.=0 //r-v,pH?2-vJ 

hp^M 

V| * * 

A reapplication of the Addition Theorem provides us a 
mathematical expression for the pressure from sphere one 
relative to sphere 2 coordinates. 

X S S ^(v 2 ’ 

p^=i i A‘ 

ri Lj Zj PiVrMzl 

f,=0 sr-ti 

hpjjid P^’(0„(^j£2^''“-(02,.«>„) 

^ V2 
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Figure 1. Orientation of sphere's origins 


C. APPLICATION OF BOUNDARY CONDITIONS 

Recalling the boundary condition for the sphere's 
surface velocity as stated in Equation (1) we obtain 



After the application of the boundary conditions, we 
use noinnal mode analysis (orthogonality) to eliminate some 
terms. In Equation (3) we note that the non-zero terms are 
those for which tj = n^, , tj = h^, and Sj = m^. 

Similarly, in Equations (4) and (5) the non-zero terms are 
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those where - n^, = m^, = n^, and JU^ = Then by 

using the relationship with sound speed, angular frequency, 
and wave niimber (c = (O/k) and dividing by a coipmon term we 
obtain 


ip^cV”' 

— - — = S S + 


h'„(ka) 


~ A ' f (ka) 
t2~0 Sz'—ti {kci) p,=^ni| 

h„{M)aTHenA) 


-ipxV”' 

° S S + 

°° A y j „ (^) / \ 

S X TT T X d P2^tiiTniiS\)' 

?.=0 “/l „^(^«)p2=?rn2l 

h„ {kd )£2r'^Gr-6'.>^+Aj 


Note that the first s\ainmation in each of these last two 
equations involves infinity. For computational purposes we 
must truncate to some finite number (JC) . Note that when K=0 
the sphere is vibrating in the monopole mode. For K=1 
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dipole modes are kept, and when K=2, quadrupole modes are 
present, and so on to higher modes. 

D. SOLVING FOR THE AMPLITUDE COEFFICIENTS 

We can write these last two equations in a matrix form 
by first constructing a large column vector whose components 
are the stacked amplitude coefficients for spheres one and 
two. The ordering of the coefficients for the spheres are 
provided below, where a"= [A„ .... ]^ is the vector of 

coefficients for sphere one. 



The C vectors above each consist of one non-zero term 
as represented by the left-hand side of Equations (7) and 
(8) . It can be shown that the terms of the and ^ 
matrices are of order (ka)^ for small ka values. This means 
that we can solve for the A-coefficients because we can 
approximate the inverse of the first large matrix of 
Equation (9). By the following development for small ka we 
can evaluate the coefficients. We drop terms on the order 
of ka to the sixth power. ■ 
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(10) 


I -B' 

B"lJ ~L-B^ I 

which leads to 

A* = C* - B'C^; - B'c^ 



Rewriting this approximation, shown below in Equation 
(11), we then have an explicit representation for the A- 
coefficients for any mode of sphere one. We can write a 
similar equation for sphere two. Note that the first term 
is due to its own mode of vibration; the next term is due to 
sphere two's monopole mode; the next three terms are due to 
sphere two's dipole modes; and the final five terms are due 
to sphere two's quadrupole modes. The Kronecker delta 
symbol {S) equals one when the n and m values equal the 
value of the second subscript (i.e. first term exists when n 
= n^ and m = m^) . In Equation (11) , I have dropped the terms 
on the order of (ka) to the sixth power. 
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E. SELF AND MUTUAL RADIATION IMPEDANCE 

The approximation for the amplitude coefficients, 
Equation (11), plays an important part of my derivation of 
the radiation impedance. We will also make use of a 
relationship developed by New and Eisler (Ref. 8). 
Radiation impedance is defined as the ratio of force exerted 
by the radiator on the median to the velocity of the 
radiator. If the velocity of the radiator is spatially 
independent, the calculation is straightfoin^rard. However, 
when the velocity does have some spatial dependence then we 
must use the equation (’ means complex conjugate) 


Z, = -^ l!p(ri)vj(ri)dSj where v, = V,P{rj) 

' V iV i 

of the 
sphere 


( 12 ). 


In applying Equation (12), note that the two identical 
spheres may be radiating in two distinct modes. Using 
Equation (12) to determine the radiation impedance for 
sphere one in the presence of sphere two and omitting terms 
on the order of (ka)* we obtain 


fe )mode « - 


* 


v:: 


‘ 'K(.kd)Sir(.6nA.) 


( 13 ). 
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Reference 9 provides an expression for the surface 
integral of the associated Legendre polynomials given in 
Equation (13). Using this expression and recalling the 
correct representation for dS in spherical coordinates, we 
obtain a general form for this surface integral 

J| ||£2:(0.«>)|fds = +1 j (14).. 

We can evaluate Equation (13) for several different 
combinations of modes (e.g. monopole-monopole, dipole- 
monopole, etc...) and determine the total radiation iitpedance 
for one sphere due to both spheres. New and Eisler 
represent the total radiation impedance of the first sphere 
into the form +Z^^*V^/V^, where is the self¬ 

radiation impedance and Z^^ is the mutual radiation 
impedance., For the monopole to monopole case we obtain the 
following expression for the mutual radiation impedance. 

Upon evaluating several combinations of modes and 
examining the calculated mutual radiation impedance I was 
able to determine an approximate form for the mutual 
radiation impedance. For~ two identical spheres each 
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radiating in a single mode of vibration the mutual radiation 
impedance is written as (I omit terms on the order of (ka)^) 



r 


!>=|n2-n,| 


(16). 


Similarly the general form for the self-radiation 
impedance for our sphere radiating in a general mode can be 
derived. Equation (17) below provides the general mode 
self-radiation impedance (I omit terms on the order of 
(^ca)^) . 



-(n, + m,)! h^(ka) 



(17) 


These two expressions provide a quick approximation for 
the mutual and self radiation impedances between two spheres 
each radiating in some general mode. I will hereafter refer 
to the mutual and self impedances calculated using Equation 
(16) and (17) as the "Modal Pritchard" approximation. Note 
that the mutual radiation impedance expression considers 
orientation between the two spheres. The following chapter 
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will investigate the utility of the mutual 
impedance approximation. Equation (16). 


radiation 




III. UTILITY OF THE MODAL PRITCHAHO APPROXIMATION 

This chapter will investigate the utility of the Modal 
Pritchard Approximation for two types of problems. The 
first problem involves looking at the interaction between 
two identical spheres. Three cases of two-body problems are 
studied: monopoles, dipoles, and quadrupoles. The second 
type of problem involves one specific three-body where we 
find a limitation of our Modal Pritchard approximation due 
to its inability to handle scattering. Throughout this 
chapter I will use a ka value equal to one. 

A. MUTUAL RADIATION IMPEDANCE 

1. Monopole to Monopole Case 

There are an infinite nximber of combinations of mode¬ 
mode interactions to consider, even for the case of two 
identical spheres. For example, considering only up to the 
quadrupole radiating modes we have 9x9=81 combinations of 
modes between two radiating spheres to examine. I will 
limit the following analysis to just three specific cases. 

The first case will be an examination of the mutual 
radiation impedance between two spheres radiating as 
monopoles, also known as the "breathing mode". The result 
obtained using the "Modal Pritchard" approximation will be 
compared to what I will call the "Simple Pritchard" 
approximation. Named after R.L. Pritchard, the Simple 


19 



Pritchard approximation was developed to provide an 
expression for the mutual radiation impedance for two 
cylindrical sources placed on an infinite baffle, and was 
further simplified by assuming a small ka value (Ref. 10). 


Z„ = -Real(Z,.)i 



(Simple Pr itchard approximation) 



Here is the mutual radiation impedance between elements 1 
and 2, and is the self radiation impedance of element 1. 
Despite the ’fact that this approximation was made for 
cylindrical piston sources on a theoretical infinite baffle, 
other researchers have used Pritchard's approximation for 
study of other shaped sources (Ref. 8) . In doing so, the 
self radiation irttpedarice is taken to be the free-field 
value. For the three-dimensional case. Equation (18) can be 
written using the spherical Hankel function as 

Z.. = Real|z,J«h.(M) (19)- 


I will compare our Modal Pritchard approximation with 
the Simple Pritchard approximation and note the differences. 
In applying Equation (19) or Equation (18), I calculated the 
self-radiation impedance by considering a sphere in free 
field environment. 

To gauge the quality of the Modal and Simple Pritchard 
approximations, I have chosen to cort^are the results of 
these methods against a strict application of the Addition 
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Theorem. Please see Appendix A for the description of my 
application of the Addition Theorem along with the computer 
codes used to generate my results. I choose to use MATLAB 
for my calculations because of the user-friendly nature of 
this engineering tool. I validated my Addition Theorem 
results with Fortran codes written by Professor Scandrett of 
the Naval Postgraduate School. Professor Scandrett has 
continually improved his code over the last eight years, 
using it for several research efforts (Ref. 11). 

Noting the appearance of infinity in Equations (4) and 
(5), we are required to truncate to some finite value for 
practical calculations. Mutual radiation impedance 
calculations using the "Addition Theorem" method have 
truncated this value to six. Even with this value, the 
computer code must solve a large (98 by 98) matrix of 
equations. 

Figure (2) shows the calculated real part of the mutual 
radiation impedance (resistance) between two monopole 
radiating spheres. Similarly, Figure (3) shows the imaginary 
part (reactance). Note first of all that the x-axis of each 
plot is the ratio of the center-to-center distance between 
the spheres to the radius of the two (identical) spheres. 
Additionally, the impedance values have been scaled by 
Ana^pc. 
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For both• the resistance and reactance parts of the 
impedance we observe that the results of the Modal Pritchard 
approximation are in excellent agreement with results of the 
Addition Theorem method throughout the plot ranges. The 
results for the Simple Pritchard approximation do not agree 
with the Addition Theorem results, especially at smaller 
values of (d/a), which corresponds to closely-spaced 
transducers. An analysis of the errors between the Modal 
Pritchard results and the Addition Theorem results, and 
between the Simple Pritchard results and the Addition 
Theorem results, was conducted. Figure (4) shows the 
results of this analysis. For ka = 1, the Simple Pritchard 
approximation breaks down when the rstio between the center 
to center distance and the radius is below 7.5. 



Figure 2. Monopole to Monopole Resistance 
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Figure 3. Monopole to monopole reactance 



2. Dipole to Dipole Case 

I conducted an analysis of the dipole-to-dipole mutual 
radiation impedance, similar to the above analysis I took 
the axes of the dipoles to be oriented along the axis 
joining the sphere centers. In terms of multipolar 
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components, this corresponds to the n=l, m=0 mode for each 
sphere. For the Simple Pritchard calculation I first 
multiplied the real part of the free-field self-radiation 
impedance for the hard sphere (oscillating in a dipole mode) 
with the appropriate Hankel fxanction, as represented in 
Equation (19). This result was then multiplied by three. 
The reason for this correction factor is that the 
conventional Simple Pritchard approximation (Equation (19)) 
does not account for any angular orientation of the two 
spheres. The factor of three correction was derived by 
comparing the Simple Pritchard Equation (19) with the full 
Addition Theorem, in a small ka limit, for the dipole 
configuration examined (Ref. 12) . For this particular 
cinalysis, I am only considering one combination (n^=l, in^=0, 
n 2 =l, and m2=0) . There are a total of nine (3x3) different 
combinations of dipole to dipole radiation cases. 

Figures (5) and (6) show the real and imaginary parts 
of the mutual radiation impedance for this dipole-to-dipole 
case. Again, note the good agreement of the Modal Pritchard 
results with the Addition Theorem results. Notice how the 
errors for the Simple Pritchard approximation increase as 
the ratio of sphere's spacing to the radius decreases. 
Figure (7) shows the errors between the results for both the 
Modal and Simple Pritchard approximations as compared to the 
results for the Spherical Addition Theorem. While the 



errors for the Simple Pritchard approximation are not too 
bad for large values of d/a, they become significant when 
d/a is less than 8.0, for Jca = 1. 
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Figure 5. Dipole to dipole resistance. 
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3. Quadrupole to Quadrupole Case 


As a final exaitple, I conducted an analysis of the 
quadrupole-to-quadrupole mutual radiation impedance. Two 
linear quadrupoles (n=2, m=0) were considered, with their 
axes aligned with the axis joining the sphere centers. For 
the simple Pritchard calculation I first multiplied the real 
part of the free-field self-radiation impedance for the hard 
sphere with the appropriate Hankel function, again using 
Equation (19). In this case, however, it can be shown that 
Equation (19) needs to be multiplied by a correction factor 
of five. Again, the reason for this correction factor, as 
in the dipole case, is that the Simple Pritchard formula 
does not account for any angular orientation of the two 
spheres. For this particulat analysis, I am only 
considering one combination (n^=2, m^=0, ^^=2, and m^^O) . 






There are a total of twenty-five different combinations of 
quadrupole to quadrupole radiation. 

Figures (8) and (9) show the real and imaginary parts 
of the mutual radiation impedance for this quadrupole-to- 
quadrupole case. We continue to note the good agreement of 
the Modal Pritchard results with , the Addition Theorem 
results. Notice again how the error in the Simple Pritchard 
results increases as the ratio between the sphere spacing 
and sphere radius decreases. Figure (10) shows the errors 
between both the Modal and Simple Pritchard results as 
compared to the results using the Spherical Addition 
Theorem. While the errors for the Simple Pritchard 
approximation are not too bad at large values of d/a, they 
become significant when d/a is less than 8.0, for ka - 1. 
This continues to show how the Simple Pritchard 
approximation breaks down as the transducers' spacing 
decreases. 
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Figure (8). Quadrupole to Quadrupole Resistance 
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ctual error 


Quadrupoles: Errors from Addition Theorenoi 



Figure (10). Quadrupoles, Errors from Addition Theorem 
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B. EXAMINATION WITH THREE A BODY PROBLEM 

To further examine our Modal Pritchard approximation 
for the mutual radiation irtpedance, we considered the 
following three-body problem (see Figure 11). 
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In this problem we will set spheres one and two to 
radiate as monopoles. These two transducers will maintain 
the same physical dimensions and operating characteristics 
as before {i.e. frequency = 474 Hz, sound speed of 1490 
m/sec, radius = 0.5 m, ka = 1, therefore the wavelength 
equals %) . Sphere three will be considered acoustically 
hard. We examine the mutual radiation iit^edance for sphere 
one due to the presence of all three spheres. Additionally, 
we will examine the effect of moving the third sphere 
horizontally away from the axis of the other radiating 
spheres. We will investigate the performance of our Modal 
Pritchard approximation for this three-body problem and 
compare these results with the results using the full 
Addition Theorem. As shown in Figure (11), sphere three's 
starting distance from the axis is one quarter of a 
wavelength (D = 71/4) . 

Figure (12) shows a plot of the mutual radiation 
resistance, comparing the results of our modal Pritchard 
approximation with those of the Addition Theorem. The x- 
axis for this plot, as well as Figures (13) and (14), is the 
ratio of the distance (P) to the radius of the spheres. 
Additionally, the results have been scaled by Ana pc. 
Figure (13) shows the reactance (or imaginary) portion of 
the mutual radiation impedance. Figure (14) shows the 
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relative error in the magnitudes of the results for the 
Modal Pritchard approximation against those for the full 
Addition Theorem. Discussions explaining my results follow 
Figure (14). 







3-Body Mutual Reactance 
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Figure (13). Mutual Reactance, Three-Body Problem 
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The modal impedance approximation as T am using it here 
neglects scattering. This is why we see the straight lines 
for the Modal Pritchard results in Figures (12) and (13). 
Because the Modal Pritchard approximation does not account 
for the scattering, this provides an argriment for using the 
full Addition Theorem to calculate impedance when scattering 
is inportant. Appendix B provides the method and corrputer 
code I used to calculate the Addition Theorem results for 
this three-body problem. 

I defined the relative error, as plotted in Figure 
(14), as (Modal Pritchard-Addition Theorem)/Addition 
Theorem. Therefore, dips in Ficfure (14) represent 
occurrences when the Addition Theorem results have a larger 
magnitude than the magnitude of the Modal Pritchard 
approximation results. Conversely, peaks represent 
occurrences for which the Addition Theorem results are 
smaller in magnitude. Through the use of the Addition 
Theorem, one can show that these peaks and valleys result 
from the effect of scattering. 

The data used to generate Figure (14) provides the D/a 
values where these valleys and peaks occur. The first two 
valleys occur at the following D/a values: 2.57 and 6.07. 
The first two peaks occur at the following D/a values: 4.57 



I propose and will sketch out how these valleys and 
peaks occur, due to scattering effects, by evaluating the 
sxim of the pressure waves arriving at sphere one, I will 
limit this evaluation by only looking at the monopole case. 
Since the Addition Theorem correctly handles scattering, the 
Addition Theorem results provide local maximum values when 
the scattering waves from sphere three are in phase with the 
direct path waves from sphere two. Similarly, the Addition 
Theorem results have local minimum values when the 
scattering waves are 180° out of phase with the direct 
path's phase. Our Modal Pritchard approximation does not 
handle the scattering effects, note the flat plots in 
Figures (12) and (13) for the Modal Pritchard results. 

Using the Addition Theorem and considering only the 
monopole case, the pressure waves are of the form of 
outgoing spherical Hankel fxinction as shown below. 

( 20 ) 

An examination of the direct path wave from sphere two 
to sphere one and of the reflected wave from sphere two, 
reflected off sphere three towards sphere one, was completed 
by Professor Scandrett (Ref. 12) . The expression for the 
sum of these waves at sphere one is provided below. The 
first term in the bracket is due to. the direct path and the 
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second term is due to the reflected (scattered) wave. The A 
is some arbitrary amplitude coefficient. 







A plot of Equation (21) is provided in Figure (15) in 
terms of D/a versus magnitude. The peaks of Figure (15) 
match the valleys of Figure (14),. Similarly, local minimums 
of Figure (15) match the peaks of Figure (14). The first 
two peaks of Figure (15) are at U/a values of 2.57 and 6.07 
and the first two local minimums are at D/a Values of 4.57 
and 7.82. The maximum values occur when the scattering 
waves are in phase with the direct path waves. The minimum 
values of Equation (21) occur when the scattering wave is 
180 degrees out of phase with the direct path. 



Figure (15): Plot of Equation (21) 
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Figure (16) provides an additional drawing of this 
three-body problem and shows the distance parameters and 
dij) used in Equation (21). 



Again, the full Spherical Addition Theorem is 
recommended over our Modal Pritchard approximation when 
scattering is an irrportant consideration as in this three- 
body problem. 
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IV. NUMERICAL MODELING WITH ATILA 

This chapter describes the ATILA finite-element model 
used to calculate the scattering of acoustic waves from an 
elastic sphere in water. The goal is to compute the so^ 
called T-matrix, which relates the scattered waves to. the 
incident waves, in an orthogonal (spherical harmonic) basis 
set (Ref. 11). 

A. MODEL INTRODUCTION 

A few comments about ATILA are appropriate at this 
time. As mentioned in the introduction, ATILA is a French 
designed finite element code specifically developed to model 
sonar transducers. ATILA, written mostly in standard 
FORTRAN 77, is the result of many years of research done at 
the Insitut Superieur d'Electronique du Nord. It is quite 
large, consisting of over 200,000 lines of code (Ref. 13). 

This tool has been used by NFS faculty and masters 
thesis students for several years and has been updated at 
least three times. Recent studies have helped to identify 
needed modifications to the code, such as increasing the 
maximiam allowable number of degrees of freedom, and 
modification of the numerical techniques used to handle the 
radiation boundaries. For example, during the calculation 
of the entries for the transition matrix (T-matrix) there 
were noted errors when compared to answers derived from a 
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known analytical solution to the problem (Ref. 11). The 
ATILA version I used (5.1.1) did increase the maximum 
allowable degrees of freedom, but has yet to address the 
radiation boundary problem. 

B. ATILA AND MESH DESCRIPTION 

1. About the Model 

A distinguishing characteristic of ATILA is the use of 
isoparametric elements for both the shell elements and the 
fluid elements. Isoparametric elements are elements that 
use the same shape function for geometric mapping as well as 
field calculations (Ref. 4). The major advantage of 
isoparametric elements is that they Tend themselves well to 
ntamerical integration. This is because of the nature of 
isoparametric elements. These elements are defined in the 
natural domain to be normalized, thus making nvimerical 
integration much easier to apply (Ref. 4). Also, when 
isoparametric elements are used, fewer nodes are needed to 
define the mesh than without such elements. 

Other nximerical techniques for solving the acoustic 
equations have been used by U.S. Navy labs. One method used 
for many years at the Naval Ocean Systems Center is called 
CHIEF (Combined Helmholtz Integral Equation Formulation) 
(Ref. 2) . CHIEF, which models acoustic radiation from 
bodies with arbitrary shapes does not use isoparametric 
elements. Thus we should e3q>ect a higher degree of accuracy 
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with ATILA over CHIEF. One advantage to CHIEF is that the 
radiation boundary condition is "built-in" to the boxandary 
element by its use of the Helmholtz Integral Equation. 
Currently, this provides an advantage over ATILA, which 
employs so-called "monopolar" and "dipolar" "dampers". 
These dampers can not represent a perfectly absorbing 
radiation boundary for a higher order multipolar field. 

The basic set of equations that allow for a large 
number of analysis to be performed include elasticity in the 
structure, the Helmholtz ecjuation in the fluid, and 
Poisson's equation in the elastic or piezoelectric material. 
Primary variables found using ATILA include the displacement 
field n in the whole structure, the electric potential 0 in 

the piezoelectric material, and the pressure P in the fluid. 
In matrix form the equations are as follows (Ref. 13). 


{[kJ-(0^\m ]) 

[is:J 

-[lI 

'll 


> ' 

[is:.,] 

[^:J 

fo] 


zz 





P 




[oF 

([h ]-<y%.]) 





( 22 ) 


Vectors in the above equations are: 

U= nodal values of the components of the displacement field, 
^ nodal values of the electric potential, 
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P= nodal values of the pressure field, 

F= nodal values of the applied forces, 

5 = nodal values of the electric charges, and 

nodal values of the integrated normal derivative of the pressure on the 
surface boundary S. 

Matrices in Equation (22) include 

= stiffiiess matrix, 

= piezoelectric matrix, 

dielectric matrix, 

M = consistent mass matrix, 

= fluid (pseudo-) stiffriess matrix, 

M\ = consistent (pseudo-) fluid mass matrix, 

L = coupling matrix at the fluid structure interface, and 

O = zero matrix. 

In addition, 

r»= Angular frequency, 

/? = Ruid density, 

c = fluid sound speed, and 

^ = means transposition of a matrix. 

2. Three-dimensional Fluid Mesh 

To run this and other finite element models you need to 
split the region under study into elements; the ATILA code 
has an extensive library of elements available. These 
elements are distinguished by a set ordering of nodes and 
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the node coordinates. Figures 16 provides a Mercator 
projection of the eight "super-elements" used for one of the 
several layers under study. These super elements are the 
start of the mesh generation. The solid lines and nodes 
(dots with numbers) define the super-elements. The dotted 
lines show how the super-elements are sub-divided. Laid out 
3-dimensnionally, Figure 17 represents the super-elements 
that make up the spherical transducer \mder study, a 
spherical shell with a radius of 0.5 m. Using a finite 
element mesh generator called MOSAIQUE, I built a complex 
mesh of the entire field using the super-elements, using 
special instructions for subdividing the elements into a 
large niomber of elements and nodes that will define each 
layer of the field. Figure 18 provides a detailed look of 
the field under study (the fluid outer boundary is at 5 
meters). 

There are several layers that make up the entire mesh 
some of which are full-layers and some are mid-layers. Each 
full-layer consists of 194 nodes of which 18 are nodes of 
the super-elements. Mid-layers have a total of 62 nodes. 
Starting from the innermost radius and working outward, the 
complete three-dimensional mesh consists of the following 
layers. First is the shell layer that consists of 194 
nodes. This is the actual spherical shell transducer that 
is being analyzed. Next, another 194 nodes are used to 
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describe the interface elements used to model the shell-to- 
water interface. This layer shares the same coordinate 
locations as the shell's nodes. Following this are ten 
fluid layers, the first two of which are 0.25 meter thick 
and the remaining are 1.0 meter thick. Each fluid layer 
consists of a full-layer and a mid-layer. Finally 
completing this three-dimensional mesh we have the radiation 
boundary elements. These last elements are used to 
prescribe monopolar-only or monopolar and dipolar radiation 
boundary conditions and are used to terminate the mesh. 



18 18 18 18 

Figure 17. The Super-Element (Macerator projection) 
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Figure 19. The complete Mesh 
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V. NDMBERICAL MODELING WITH ATZLA: RESULTS AND DISCUSSION 

A. OBJECTIVES AND THE T-MATRIX 

My objectives for the numerical modeling portion of 
this thesis were to investigate the finite element mesh 
mentioned earlier, simplify the process of calculating the 
Transition Matrix (or T-Matrix), and provide a new method of 
evaluating the off-diagonal elements of the T-matrix. The 
finite element code is used to numerically calculate the 
scattered pressures from an object under study. We then use 
these results to calculate the object's specific T-Matrix. 

In this chapter, I will document the results of the T- 
matrix calculation and compare with previous numerical 
modeling efforts. Several references mention and thoroughly 
discuss the definition, usage, and calculation of the T- 
Matrix (Ref. 3 and 5) . I will simply mention here in the 
text that the T-Matrix describes the scattering properties 
of a specific body. It does this by using a discrete basis 
set spherical harmonics in our case, and is used to show how 
the incident pressure waves translate to scattered waves 
from the body. Please see Appendix C for more information 
about the T-Matrix and how I calculated the elements of the 
T-Matrix for our radiating body. 
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B. RESULTS 

In order to gauge the accuracy of our T-Matrix 
calculation we compare our answers with an, analytically 
derived equation. An expanded discussion of the derivation 
of the analytical solution can be found in a Master's Thesis 
by Ruiz (Ref. 5 and references therein). For a thin 
spherical shell the T-Matrix is diagonal, and are given by 


iZ.jSka) + p,c,jSka) 
iZ.hAka) + P,C,hAka) 


( 23 ) 


where 


,Ac,A f[Q^-2(l+v)](Q^+l-v-A.)-/l.(l+v)' 


Z, = ! 


= c 


+1 + v—A 




- 9 /' - f / Vf#!./ 1^9 I 9 '•'I**' * -fc/ 

C, ' 'Va(i-»') 

V = Poisson^s ratio (0.33); h = shell thickness (1 cm) 

Cp - shell's plate velocity 
E = Young's modulus {0.\25x\0^^) 
a = spherical shell radius (0.5 meter) 

- shell's density (7500 kglm^)\ f = frequency (474 Hz) 
p^ and Cf - fluid density (1000 kgfm^) and speed (1490 m!s). 


The calculation for this analytical solution was 
performed using MATLAB and the computer program can be found 
at the end of Appendix C. In the following discussion, when 
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referring to errors, I am comparing the calculated values 
with the analytical solutions from Equation (23). 

Using the thin spherical shell parameters listed above 
and the MATLAB computer script provided in Appendix C, I 
calculated the analytical T-Matrix diagonal elements (the 
off-diagonal elements are zero for a spherically symmetric 
scatter). These diagonal elements are provided in Table I 
below. A value of ka = 1 was used (frequency = 474 Hz). 


ANALYTICAL DIAGONAL T-MATRIX ENTRIES 

ELEMENT 

REAL 

IMAGINARY 

MAGNITUDE 

PHASE 


PART 

PART 


(DEGREES) 

Til =Roo 

-1.3465E-02 

1.1525E-01 

1.1604E-01 

96.6635 

T 22 = Ri-i 

-6.0255E-03 

7.7390E-02 

7.7624E-02 

94.4520 

T 33 = Rio 

-6.0255E-03 

7.7390E-02 

7.7624E-02 

94.4520 

T44 = Rii 

-6.0255E-03 

7.7390E-02 

7.7624E-02 

94.4520 

T 55 = R 2-2 

-6.1773E-04 

-2.4847E-02 

2.4854E-02 

-91.4242 

T66 = R 2 -I 

-6.1773E-04 

-2.4847E-02 

2.4854E-02 

-91.4242 

T 77 = R 20 

-6.1773E-04 

-2.4847E-02 

2.4854E-02 

-91.4242 

Tgs = R 2 I 

-6.1773E-04 

-2.4847E-02 

2.4854E-02 

-91.4242 

T 99 = R 22 

-6.1773E-04 

-2.4847E-02 

2.4854E-02 

-91.4242 


Table I. Analytical Diagonal T-Matrix Entries, 


I then used the finite element code (ATILA) to model 
the scattering pressure field generated by the thin 
spherical shell. With the scattered pressure results, I 
calculated the T-Matrix entries colxomn by coliamn. These 
results are provided below in Table II. 
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Numerical Diagonal T-Matrix Elements 

Element 

Real 

Imaginary 

Magnitude 

Phase 


Part 

Part 


(Degrees) 

T„ 

-1.6107E-02 

1.2894E-01 

1.2994E-01 

97.1202 

T22 

-8.0999E-03 

7.9943E-02 

8.0353E-02 

95.7855 

T33 

-7.7968E-03 

7.9777E-02 

8.0157E-02 

95.5819 

T 44 

-8.0801E-03 

7.9946E-02 

8.0353E-02 

95.7713 

T 55 

-2.1393E-03 

-1.5865E-02 

1.6009E-02 

-97.6794 


-2.2524E-03 

-1.5923E-02 

1.6082E-02 

-98.0513 

T77 

-2.1536E-03 

-1.5930E-02 

1.6075E-02 

-97.6994 

Tgs 

-2.2501E-03 

-1.5925E-02 

1.6084E-02 

-98.0420 

T99 

-2.1498E-03 

-1.5860E-02 

1.6005E-02 

-97.7193 


Table II. Nxamerical Diagonal T-Matrix Elements. 


Table III, next page, provides an error analysis for 
the current results. Unfortunately, the monopole result 
(i.e. the first diagonal element) is much worse than 
previous results. All other diagonal elements show some 
improvement over previous work, between 6 and 56 percent 
improvement in their magnitudes. But the most important 
indicator is the monopole element; therefore, our new mesh 
does not provide a significant improvement to the T-Matrix 
calculation. As the ATILA code is further improved we hope 
that these errors will be corrected. 
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Error Comparison Current vs. Previous 

Element 

Magnitude 

Phase 

Previous 

Previous 

Magnitude 

Phase 


Error 

Error 

Mag. Error 

Phase Error 

Improve 



(percent) 

(Degrees) 


(Degrees) 

(percent) 

(Degrees) 

Tn 

11.984 

0.457 


-0.360 

-1216.94 

- 0.10 

T22 

3.515 

1.333 

6.490 

0.329 

45.84 

- 1.00 

T33 

3.263 

1.130 

7.330 

-0.500 

55.48 

-0.63 

T44 

3.516 

1.319 

6.470 

0.314 

45.66 

-1.01 

T55 

-35.588 

-6.255 

-37.800 

-9.186 

5.85 

2.93 


-35.296 

-6.627 

-38.000 

-9.006 

7.12 

2.38 

T77 

-35.324 

-6.275 

-37.700 

-7.823 

6.30 

1.55 

Tgs 

-35.288 

-6.618 


-8.986 

7.14 

2.37 

T 99 

-35.605 

-6.295 

-37.900 

-9.106 

6.05 

2.81 


Table III. Error Comparison Current vs. Previous. 

The T-Matrix must be diagonal due to the symmetry of 
the object. That means that the off-diagonal elements must 
be zero. Previous work dociimented eight significant off- 
diagonal elements. Current results reveal these same eight 
elements are the most significant off-diagonals. However, 
the new results show improvement, meaning the new off- 
diagonal values are much closer to zero than previous 
calculations. Along with matching the analytical diagonal 
elements, making all off-diagonal elements closer to zero is 
another indication of how well our finite element method is 
working. Table IV offers a comparison of this improvement. 
The last coliimn of this table shows the percentage closer to 
zero the new off-diagonal magnitudes are compared to the 
previous work. 
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SIGNMCANT OFF-DIAGONALS 


Element 

Real 

Part 

Lnaginary 

Part 

Improvement 
(percent closer to zero) 

T 71 

3.21650E-04 

6.87400E-04 

18.96 


-9.89950E-08 

-1.78470E-07 

99.12 

T84 

-3.02040E-18 

8.25020E-19 

100.00 

T 95 

1.18910E-07 

-9.80070E-07 

95.46 

T 26 

3.69920E-06 

1.98190E-06 

62.58 

Ti7 

-5.87890E-04 

1.30000E-03 

79.74 

T 48 

9.56520E-07 

2.56450E-07 

98.07 

T 59 

2.68580E-04 

1.18460E-(H 

97.69 


Table IV. Significant Off-Diagonals. 


G. FURTHER AIU^YSIS AND DISCUSSION 

Some Matrix algebra is employed to further characterize 
how well the T-Matrix off-diagonals were calculated. The 
following equation provides a normalization method used to 
characterize the off-diagonal elements relation to one. 


1—- 
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T 
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0 - 0 

V'^99 


( 24 ) 


By multiplying both sides of the T-Matrix with the 
diagonal matrix as described in Equation (24) I will obtain 
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a matrix with one's along the diagonal and relative errors 
in the off-diagonal elements. By examining the magnitudes 
of the most significant off-diagonal elements, I can quickly 
gauge how "well" I calculated the T-Matrix. Again we see 
the same eight off-diagonals show the most significant 
levels with respect to the normalization. Although they are 
the most significant they are small. This method of 
analyzing the T-Matrix's off-diagonals provides a "quick- 
check" as to how well our niomerical analysis is working. 
Table V provides our results for this analysis. Notice how 
small the magnitudes of these elements. All other elements 
have magnitudes on the order of element T,^ and smaller. 


T-MAHUXNORM/ 
(MOST SIGNMCA 

OJZATION ANALYSIS 

NT OFF-DIAGONALS) 

Element 

Real 

Imaginary 

Magnitude 


Part 

Part 


T71 

6.96163E-03 

1.50759E-02 

1.6606E-02 

T 62 

-2.65519E-06 

-5.0I816E-06 

5.6773E-06 

T 84 

-8.44548E-17 

2.12802E-17 

8.7095E-17 

T 95 

5.96808E-05 

1.55646E-05 

6.1677E-05 

T26 

1.01798E-(M 

5.71578E-05 

1.1675E-(M 

Ti7 

-1.30097E-02 

2.89475E-02 

3.1737E-02 

T48 

2.64606E-05 

7.65933E-06 

2.7547E-05 

T59 

-9.58209E-03 

1.56365E-02 

1.8339E-02 


Table V. T-Matrix Normalization Analysis. 

I 
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VI. CONCLUSION AND FOLLOW-ON PROPOSALS 

A. CONCLUSION 

I , 

I have derived an approximation for the self and mutual 
radiation impedances for spherical transducers of small ka 
values by considering two spheres in isolation. In this 
thesis, I have decided to call this approximation "Modal 
Pritchard". This approximation was tested for a general 
sphere of radius a=0.5 meters, frequency of 474 Hz, in water 
with sound speed of 1490 m/sec. Two general problems were 
tested; (1) two bodies radiating in three modes (monopoles, 
aligned dipoles, and aligned-linear quadrupoles) and (2) a 
three-body problem in which two bodies radiated in an array 
and the third acts as a moving scatterer. In all the two- 
body problems, the Modal Pritchard approximation provides 
mutual radiation impedance results close, if not exactly, 
matching the full Addition Theorem results. The same cannot 
be said for the three-body problem where scattering is 
important. Therefore, when scattering is important we 
recoitimend using the full Addition Theorem. In some cases, 
our Modal Pritchard approximation may be a useful tool in 
the design phase of low-frequency active sonars. 

The nxamerical modeling portions of this thesis have 
updated and streamlined the process of calculating the 
scattering characteristics and T-matrix for an object under 
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study. Using a more refined mesh, only a partial 
improvement has been realized for the T-matrix calculation 
for a thin spherical-shell radiator. Unfortiinately, current 
work showed worse errors for the monopole element of the T- 
Matrix. 

We also noted some improvements in the calculation of 
the off-diagonal elements.. Most off-diagonal elements 
showed an improvement of between 63 to 100 percent closer to 
zero as compared to previous work. I also provided a quick 
and easy means for determining how well we calculated the T- 
matrix off-diagonal elements by "normalizing" the matrix. 
With this process, we can quickly gauge the off-diagonal 
elements. AS a result, seven off-diagonal elements were 
between 0.02 and 0.000006 as compared to the number one. 
All remaining off-diagonal elements (65) were smaller than 
9.0*10'^’ (most significantly smaller). Current work shows 
our T-matrix calculation was worse than previous work for 
the monopole element, but showed some improvements for other 
elements of the T-matrix. There is still room for further 
inprovement which may be realized when ATILA radiation 
boundary problems are addressed. 



B. FOLLOW-ON PROPOSALS 


Further work with the "Modal Approximation" to the 
mutual radiation impedance include studying other problems 
such as arrays of multiple transducers and calculating the 
pressure field generated from their radiation. This 
investigation does not need to be limited to comparisons 
with computer generated results from the Addition Theorem. 
Follow-on work could use this Modal Pritchard approximation 
to compare with laboratory results of real transducers. 
With meticulous data gathering, one could study how well 
this approximation does in the "real" world. 

As mentioned earlier, there continues to reside 
problems with the radiation boundary treatment in ATILA. 
This ^problem was to have been addressed with an even newer 
version of ATILA to include a tool called EQI. In reference 
to the second portion of this thesis, it is hoped that this 
modification will further improve the T-matrix calculation. 
Currently, about half of our T-matrix diagonal elements show 
magnitude errors below 7.3 percent and the rest are about 38 
percent in error. Current phase errors are below 9 percent. 
Further analysis using the modified ATILA with EQI may reach 
our goal of relative errors below one percent. Future work 
could also use other refined meshes that may also provide 
further improvements. 
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APPENDIX A: ADDITION THEOREM CALCULATION OF MUTUAL RADIATION 

IMPEDANCE 

Here I will explain how I used the addition theorem to 
calculate the mutual radiation impedance for two identically 
hard spherical shells. I will refer to some equations in 
the text of this thesis throughout this explanation and will 
provide my MATLAB program along with some user defined 
ftinctions. 

Recall the application of the Addition Theorem to write 
the pressure from one sphere relative to the coordinate 
system of the other sphere. Equations (4) and (5) of the 

text. The application of the boundary conditions as written 

in Equation (6) brings us to Equations (8) and (9). I will 
rewrite the last two equations here. 

°° A r, J ' I \ 

E L AL , \ E a Ui’P,’?2’mi’52r 

hpjkd 

~ A 1 7 ’ / A 

E E -- E a \n 2 ’P 2 ’ti^m 2 ^Sir 

PjH'SrmJ 

hpikd 
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These two equations are put into matrix form as given 
by Equation (10). Notice the infinity sign in the 
expression above. For practical purposes we have to 
truncate this infinity to some finite number K. Using K = 
6, the left hand side of Equation (10) becomes a large 98 by 
98 matrix multiplying a 98 element column vector for the 
left hand side and the right hand side is 98-element colvimn 
vector with only two non-zero entries. We solve for the 
anplitudes in MATLAB by using a forward divide as follows. 



With K = 6 the matrices have the following sizes, 



Now we have solved the unknown A-coefficients and are 
ready to use these results to find the mutual radiation 
iitpedance. 

Using Equation (13) with our solved A-coefficient 
values, we can easily solve for the radiation impedance. 
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given that the total radiation impedance can be written as 
+Z^*VJV^. I first solve for the total radiation 
impedance by setting both V's equal to one. Then by setting 
^2 = 0 I obtain the self-radiation impedance (Z^j) result. I 
then subtract Z^^ from to obtain the mutual-radiation 

impedance (Z^j) result. 

Here are the MATLAB computer scripts used to generate 
these results along with all supporting user defined 
functions written for this analysis. 


MAIN MATLAB PROGRAM 
% Joseph L. Day 

% Mutual iitpedance using the Addition Theorem 
% 2 spheres along z-axis, with top sphere moving away 
% Last updated July 27, 1999 
format long 

K=6; % Truncated integer vice infinity 

K2=(K+1)'^2; % Useful number, it is the size of many vectors 

k=2; % Acoustic wave niomber. 

radius=0.5;ka=k*radius; % Both spheres with same radius. 
theta=0;phi=pi; % Angle from 2's origin to I's 

nl=0; ml=0; n2=0; m2=0; % Defines the mode for each of the 

two spheres 

positionl=l; % Positionl and Position2 are integers change 
position2=l; % as follows: nl=0,ml=0 then positionl=l; 

% nl=l,ml=-l then positionl=2; and 
% nl=l,ml=0 then positionl=3,..and so on. 
V1=1;V2=0; % Spheres velocity amplitudes 

% Both equal to one gives Z(total) 

% With Vl=l and V2=0 gives Zll 
speed=1490; % speed of sound in the water. 

rho=1000; % density of the water. 

templ=eye(K2);tenp2=zeros(K2,K2); 

Big=[tempi temp2;temp2 tempi]; % Matrix allocation 
count=l; % Initialize a cotinting process 

for d=2.5:.25:24 % Distance with 0.25 meter increment 
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countl=l;count2=l; % Counting process. 

Btertpld=zeros (1,K2) ; 

Btemp2d=zeros(1,K2); 

for t=0 : 6 

for s=-t:t 

ptempla=abs(t-nl) ;ptemplb=abs(s-inl) ; 
ptemplc=inax(ptempla,ptemplb) ; 
ptemp2a=abs(t-n2);ptemp2b=abs(s-in2); 
ptemp2c=inax (pteinp2a, ptemp2b) ; 
count3=l;count4=l; 
for pl=ptemplc:2:t+nl 
M=s-ml; 
if pl<0 

p=(-1)^abs(M)*fact(pl-abs(M))7 
fact(pl+abs(M))*legendre(pi,cos(theta)); 

,'else ■ 

P=legendre(pi,cos(theta)); 

end 

Btempla=a(nl,pl,t,inl,s)*shank2(pl,d); 
Btemplb=P(l+abs(M)) *exp(i*M*phi) ; 

Btemplc(1,counts)=Btempla*Btemplb; 
count3=count3+1; 
end ' 

Btempldd, countl) =s\im(Btemplc) *jprime (nl, ka) /hprime (nl,ka) 
coxantl=countl+l ; 
for p2=ptemp2c:2:t+n2 
M=s-m2; 
if p2<0 

P= (-1)'^abs (M) *fact (p2-abs (M)) / 
fact(p2+abs(M))*legendre(p2,cos(theta-pi)) 

else 

P=legendre(p2,cos(theta-pi)); 

end 

Btemp2a=a (n2, p2, t,m2, s) *shank2 (p2, d) ; 
Btemp2b=P(l+abS(M))*exp(i*M*(phi+pi )); 

Btemp2c (1, co\int4) =Btemp2a*Btemp2b; 
count4=count4+l; 

end 

Bteirp2d (1, count2) =sum (Btemp2c) * jprime (n2, ka) /hprime (n2, ka) 
count2=count2+l; 
end ■ 

end .. . ■ 
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Big(positionl,K2+l:2*K2)=Btempld; %Placing values 
Big(K2+position2,l:K2)=Btemp2d; %in the proper row. 
C(2*K2,l)=zeros; 

C(positionl,l)=-i*rho*speed*Vl/hpriine(nl,ka); 

C(position2+K2,1)=-i*rho*speed*V2/hprime(n2,ka); 
A=Big\C; 


Temporary=STain(conj (A(K2+1:2*K2,1)) ' . *Bteinpld*sbessj (nl,ka)) * 
hprime(nl,ka)/jprime(nl,ka) ; 

Zrl=4*pi*radius''2/Vl* (A{1,1) *shank2 {nl,ka) + Temporary) ; 
Zrvector(count)=Zrl; 
count=count+l; 

end 

Zrlscaled=Zrvector/{4*pi*radius^2*rho*speed); 

% This last value is the final answer and has been 
% scaled for plotting purposes. 

% Also for plotting purposes the final answer is calculated 
% for ranges (d/a) from 2.5 to 24 in increments of 0.25. 

% All plots were done using Microsoft Excel, since that 
% program works better with Microsoft Word for cutting and 
% pasting graphs. 
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Supporting user defined functions : 
function A=a (s, t, r,u,in) 

% computes the function a{s,t,r,u,m) of the Addition Theorem 
% as written in King and Van Buren 

% Written by Joseph L. Day 
numl={2*s+l)*(2*t+l)*fact(s-u)*fact(t- 
m+u)*fact(r+m)*fact((s+t+r)/2); 
deml=fact{(r+t-s)/2)*fact((r+s-t)/2)* 
fact{(s+t-r)/2)*fact(s+t+r+l); 
terml=numl/deml; * 

wmin=0.5*max([r-s-t,s-r-t-2*m+2*u,t-s-r+2*u]); 
wmax=0.5*min([s+t-r,r+t-s-2*m+2*u,r+s-t+2*u]); 
s\im=0; 

for w=wmin:wniax 

sum=sum+(-1)^w*bc(s+t-r,(s+t-r)/2+w)* 

be (t+r-s,(t+r-s)/2+m-u+w)*bc(s+r-t,(s+r-t)/2-u+w); 

end 

A=real (i^'(s+t-r) ) *terml*sum; % Output from this function 


function y=bc(n,m) 

% Computes the binomial coefficient 

% Written by Joseph L. Day 
if m<0 
y=0; 

elseif n-m<0 
y=0; 
else 

y=prod(l:n)/(prod(l:m)*prod(l:n-m)); 

end 


fianction y=fact (n) 

% Computes the factorial of n 
% Written by Joseph L. Day 
y=prod(l:n); 
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function hn=shank2(n,x) 

% Computes the Spherical Hankel f\anction of the second kind 

% Written by Joseph L. Day 

hn=sqrt(pi/(2*x))*besselh(n+.5,2,x); 


function jp=jprime(n,x) 

% Computes the first derivative of the Spherical Bessel 
% fiinction. 

% Written by Joseph L. Day 
jp=sqrt(pi/{2*x))*(besselj(n-.5,x)- 
{n+l) *besselj (n+.5,x) / (x) ) ; 


function hp=hprime(n,x) 

% Computes the first derivative of the Spherical Bessel 
% function. 

% Written by Joseph L. Day 
hp=sqrt(pi/(2*x))*(besselh(n-.5,2,x)- 
(n+l) *besselh{n+.5,2,x) / (x) ) ; 


function jn=sbessj(n,x) 

% Computes the Spherical Bessel function 

% Written by Joseph L. Day 
jn=sqrt(pi/(2*x))*besselj(n+.5,x); 
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APPENDIX B: THREE-BODY ADDITION THEOREM CALCULATION 

The solution for the mutual radiation impedance of the 
three-body problem as described in the Chapter III part B 
using the Addition Theorem is provided in this appendix. We 
start by writing the three radial velocities as expressed in 
Equation (2). The pressure from each of the three spheres 
(before we set sphere three as acoustically hard) can be 
expressed in the form of Equation (3). 

As before we utilize the Addition Theorem to write the 
pressures from spheres 2 and 3 relative to sphere one (see 
Equation (4)) . This gives us the expression for Pj' and Pj^. 
The boundary conditions on the surface of each of the 
spheres leads to the following three equations. 


v.=— 

cop ^ or rr« 

V. = —^[P2 + P"+Pl] 

(Op^ar 

^2 = —|-[P, + /^ + Pl] 

COp^ar r2=° 

Then we obtain three big equations that are similar to 
Equations (7) and (8) each with an additional term. For 


example looking at Equation (7), the additional term will be 
much like the term in the parenthesis but will involve the 
coefficients for the third sphere (A^) . 
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Following a similar approach as used in Appendix A, we 
put our Equation (7)-type expressions in matrix foimi. We 
can write the following matrix relation. 


7 

B\ B\~ 

- 






I Bl 

-> 

A 

— 


where 



Bl I _ 

- » 


.cl 


--1 

CO 


As before we will truncate the infinity value to K = 6 
gives expression of the following sizes. 


- 

"A’5 


c' 


~A^s' 




& 

147jc147 

147jc1 

— 

& 

c \ 

solving 



147jc147 

/ 

c \ 


We have solved for the three sphere coefficients and 
are now ready to complete our evaluation of the mutual 
radiation impedance. Using Equation (13) with these A- 
coefficients I can solve for the radiation impedance. Now I 
can write sphere one's total radiation irrpedance as 
+ Now I make sphere three acoustically 

hard by setting the radial velocity amplitude in my cortputer 
code to zero. Next I solve for the total radiation 
impedance by setting sphere one and two's radial velocity 
anplitudes equal to one. I store this result and then 

calculate sphere one's self-radiation irtpedance by setting 
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sphere two's radial velocity amplitude to zero. I subtract 
this last result from the total radiation iiipedance, and 
that provides my mutual radiation impedance for sphere one 
due to all three spheres. The plotted data is calculated by 
incrementing the range of sphere three as it moves away from 
the two-sphere axis. 

Here is the MATLAB computer code used to generate these 
results. All supporting functions have already been 
provided. 

MATLAB Program for 3-body mutual radiation impedance 


% Written by Joseph L. Day 

% Mutual impedance using the Addition Theorem 
% Two spheres radiating, third (hard) sphere moving away 
% Last updated August 16, 1999 
format long 


K=6; 


% Truncated integer vice infinity 


K2=(K+1)^2; 
k=2 ; 

radius=0.5;ka=k*radius; 
theta=0;phi=pi; 
nl=0; ml=0; n2=0; m2=0; 


% A useful number 
% waveNxiniber. 

% for all three spheres 
% Angle from 2's origin to I's 
% Defines the mode for each sphere 


n3=0; m3=0; 

positionl=l; 

position2=l; 

position3=l; 

V1=1;V2=1;V3=0; 

speed=1490; 

rho=1000; 


% Monopole position for matrices 
% another monopole; 

% like a monopole 
% spheres velocity amplitudes 
% speed of sound in the water 
% density of the water 


templ=eye(K2);temp2=2eros(K2,K2); 

Big=[tempi temp2 temp2;temp2 tempi temp2;temp2 temp2 tempi]; 

% setup some space for future calculation 

count=l; % initialize a counting 

process 

phicheck=0; 
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for d=pi/2:.25:24 % Distance with 0.25 meter 

increment 
dl2=pi; 

countla=l;count2a=l;count3a=l;% initialize counts. 
countlb=l;count2b=l;count3b=l; 
Btertplda=zeros(l,K2);Btempldb=zeros(l,K2); 
Btemp2da=zeros{l,K2);Btemp2db=zeros(l,K2); 
Btemp3da=zeros(l,K2) ;Btemp3db=zeros(l,K2) ; 
for t=0:6 

for s=-t:t 

ptempla=abs(t-nl);ptemplb= 

abs(s-ml);ptemplc=max(ptempla,ptemplb); 
ptemp2a=abs(t-n2);ptemp2b= 

abs(s-m2);ptemp2c=max(ptemp2a,ptemp2b); 
ptemp3a=abs{t-n3);ptemp3b= 

abs{s-m3);ptemp3c=max(ptemp3a,ptemp3b); 

count3=l; count4=l; count5=l; count6=l; cotint7=l; coxmt8=l ; 
for pl=ptemplc:2:t+nl %for B(lto2) 

M=s-ml; 
if pl<0 

P= (-1)''abs (M) *fact (pl-abs (M)) / 

fact(pl+abs(M))*legendre(pi,cos(theta)); 

else 

P=legendre(pi,cos(theta)); 

end 

Btempla=a(nl,pl,t,ml,s)*shank2(pl,dl2); 

Btemplb=P(l+abs(M))*e5{p(i*M^phi); 

Btemplc(1,count3)=Btempla*Btemplb; 

coun t 3 =coun t 3 + 1 ; 

end 

Btemplda(l,countla)=s\im(Btemplc) *jprime(nl,ka) / 

hprime(nl,ka); 

countla=countla+l; 
phil3=pi/2+atan(pi/2/(d)); 
thetal3=3*pi/2; 

for pl=ptemplc:2:t+nl %for B(lto3) 

M=s-ml; 
if pl<0 

P=(-1)^abs(M)*fact(pl-abs(M))/ 

fact(pl+abs(M))*legendre(pi,cos(thetal3)) 

else 

P=legendre(pi,cos(thetal3)); 

, • end 

dl3=sqrt((pi/2)^2+(d)^2); 

Btempla=a(nl,pl, t,ml, s) *shank2 (pl,dl3) ; 
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Btemplb=P(l+abs(M) )*exp(i*M*phil3) ; 

Btemplc (1, count4) =Btertpla*BteInplb; 
count4=count4+l; 

end 

Btempldb (1, count lb) =stun(Btemplc) * j prime (nl,ka) / 

hprime(nl,ka); 

countlb=countlb+l; 

for p2=ptemp2c:2:t+n2 %for B(2tol) 

M=s-m2; 
if p2<0 

P= (-1) ^abs (M) *fact (p2-bs (M)) /fact (p2+abs (M)) * 
legendre{p2,cos(theta-pi)); 

else 

P=legendre(p2,cos(theta-pi)); 

end 

Btemp2a=a(n2,p2,t,m2,s)*shank2(p2,dl2); 

Btemp2b=P(1+abs(M))*exp (i*M*(phi+pi)); 

Btemp2c {1, co\ant5) =BteIrp2a*Btemp2b; 
count5=count5+l; 

end 

Btertp2da(1,count2a) =sum{Btemp2c) *jprime(n2,ka) / 

hprime(n2,ka); 

count2a=co\mt2a+l; 
phi23=pi/2-atan(pi/2/(d)); 
theta23=3*pi/2; 

for p3=ptemp3c:2:t+n3 %for B(2,3) 

M=s-m3; 
if p3<0 

P= (-1)'^abs (M) * f act (p3-abs (M)) / 

fact(p3+abs(M))*legendre(p3,cos(theta23)); 

else 

P=legendre(p3,cos(theta23)); 

end 

d23=sqrt( (pi/2)''2+(d)^2) ; 

Btemp3a=a (n3 ,p3, t ,m3, s) *shank2 (p3, d23) ; 

Btemp3b=P {1+abs (M)) *e3{p (i*M* {phi23)) ; 

Btemp3c (1, counts) =Bterrp3a*Btemp3b; 
count6=count6+l; 

end 

Btemp2db(l,co\int2b)=s\am(Btemp3c) *jprime(n3,ka) / 

hprime {n3,ka) ; 

count2b=count2b+l; 

for p3=ptemp3c:2:t+n3 %for B{3,1) 

M=s-m3; 
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if p3<0 

P=(-1)^abs(M)*fact(p3-abs(M))/ 
fact(p3+abs(M))*legendre(p3,cos{thetal3-pi)) 

else . ' ' '. 

P=legendre(p3,c6s(thetal3-pi)); 

end 

d23=sqrt((pi/2)''2+(d)''2); 

Btemp3a=a(n3,p3,t,in3,s)*shank2(p3,dl3); 

Btemp3b=P(l+abs(M))*exp(i*M*(phil3+pi)); 

Btemp3c (1, count?) =Bteitp3a*Bteinp3b; 
count7=count7+1; 

end 

Btemp3da (1, count3a) =siam(Btemp3c) * jprime {n3, ka) / 

hprime{n3,ka); 
count3a=count3a+l; 
for p3=ptemp3c:2:t+n3 %for B(3,2) 

M=s-m3; 
if p3<0 

P=(-1)^abs(M)*fact(p3-abs(M))/ 
fact(p3+abs(M))*legendre(p3,cos{theta23-pi)) 
else 

P=legendre(p3/COs(theta23-pi)); 

end 

Btemp3a=a{n3,p3,t,m3,s)*shank2(p3,d23); 
Btemp3b=P(l+abs(M) )*exp(i*M*(phi23+pi)) ; 

BteInp3c (1, coTontS) =Bterr^3a*Bteinp3b; 
count8=count8+l; 
end ^ 

Btenp3db(l,coiant3b)=sum(Bteinp3c) *jprime(n3,ka) / 

hprime(n3,ka);. 

count3b=count3b+l; 

end 

phicheck=phicheck+.25; 
end 

Big(positionl,2*K2+l:3*K2)=Btenrpldb; % Placement • 

Big(positionl,K2+1:2*K2)=Btemplda; 

Big(K2+position2,l:K2)=Btemp2da; % of vectors. 

Big(K2+position2,2*K2+1:3*K2)=Btemp2db; 
Big(2*K2+position3,1:K2)=Btemp3da; 

Big (2*K2+position3 , K2+1 : 2*K2 ) =Btemp3db; 

C(3*K2,l)=zeros; 

C (positionl,1)=-i*rho*speed*Vl/hprime(nl,ka); 
C(position2+K2,l)=-i*rho*speed*V2/hprime(n2,ka); 

C(position3+2*K2,1)=-i*rho*speed*V3/hprime(n3,ka); 
A=Big\C; 



Temporary=suin{conj (A(K2+1:2*K2,1)) '.*Btemplda*sbessj {nl,ka)) 
*hprime(nl,ka)/jprime(nl,ka) ; 

Teinporary2=STam(conj (A(2*K2+1:3*K2,1)) '.*Btempldb*sbessj (nl,k 
a))*hprime(nl,ka)/jprime(nl,ka); 

Zrl=4*pi*radius^2/Vl*(A(l,1)*shank2(nl,ka)+Temporary+Tempora 
ry2) ; 

Zrvector (count) =Zrl; 
cotint=count+l; 

end 

Zrlscaled=Zrvector/{4*pi*radius''2*rho*speed) ; 


73 



(THIS PAGE INTENTIONALLY LEFT BLANK) 



APPENDIX C: T-HATRIX CAIiCULATION 

Upon obtaining the scattered pressures from the ATILA 
code for the user defined mesh, I then calculated the 
Transition Matrix (T-Matrix) for the user-defined 
transducer. For this analysis we used a thin spherical 
shell for our transducer. 

Before I spell out my method for calculating the T- 
matrix, I should provide the reader with some background 
concerning this T-matrix (Ref. 3) . First suppose the 
incident pressure to an object are standing spherical waves 
as represented by 

We use this T-matrix to transform from our incident 
pressure to the scattered pressure from the object subjected 
to this incident pressure. Basically, the T-Matrix 
describes the scattering property of each unique scatter (or 
in this case transducer) by using a discrete basis set of 
spherical harmonics. I determine the T-Matrix for my 
spherical shell by analyzing the scattering characteristics 
of this shell subject to a "standing" wave. . Mathematically 
speaking, the T-matrix transforms our incident standing 
waves to scattered outgoing waves in the form of Hankel 
functions shown below. 
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p{r,e,(j>)=}ifir)^^{e,(p). 

The following equation indicates the utility of the T- 
Matrix. The values represent the reflection coefficients 
from the spherical shell and constitute elements of the T- 
Matrix. 

T-Matrix / v 

n ^ . , « , 

'---,-' ... '-—V-' 

S tsn. ding Wave Radiating Wave 

Here is how I used these ideas to develop my T-matrix 
calculation. Starting with the output data from ATILA which 
is basically N-nxmiber of data points in the following form: 
for i = 1 to N we have x^, Yi> and P^. Then for all N- 

points I need to calculate the following product 

I will do this product up to the quadrUpole, that is 
nine times for each n and m combination (n=0 and m=0, n=l 
and m=-l, ... n=2 and m=2) . So for example if we have 100 

data points, the number of calculated complex numbers is 900 
(100 data points up to the quadrupole) . The resulting 
matrix is of the form: 
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calculated complex numbers 

■ t ■ 

A 

{unknowns) 


■ t ' 

Vxessure 

(fromATILA) 

= 

100 9 

(9^1) 

i 


(100a:1) 

i 






We solve for the iinknovm A-coefficients in MATLAB by 
the following equation [A] = [calculated numbers]\[P]. The 
forward divide sign in MATLAB is a solution for the A vector 
in a least squares sense to this over-determined system of 
equations [calculated n'umbers] * [A] = [P] . The effective 
rank of the calculated niimbers matrix is determined from the 
QR decomposition with rank pivoting. The relationship of 
the T-matrix with the incident pressure amplitude vector [B] 
and the A-coefficients is as follows 



Boo 


Aoo 



Bo-i 


Ao-i 

T-matrix 


Bio 

— 

Aio 



B22 


■ ^ 

_1 
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Since we obtain bur results from ATILA for each n and m 
combination we solve for each colximn of the T-matrix in the 
following way. For exanqple, we obtain ATILA's first 
scattering pressure results when we model the incident 
pressure with the standing waves using n=0 and m=0. We then 
use the scattering pressures to solve for the A- 
coefficients. Then we solve for the first column of the T- 
matrix by solving for the following equation. 


1 ■ 


Aoo 

0 


-Ao-i 

0 

= 

Aio 

0 


_i422 _ 


Proceeding this way we can solve for all 9-columns of 
the T-matrix by changing the n and m values Of the incident 
pressure accordingly. One final comment about this 
procedure to calculate the T-Matrix concerns the form of the 
incident wave. Because of certain limitations with the 
ATILA code we used a "traveling" wave vice the "standing" 
wave. The form of the incident pressure wave possible is a 
Hankel function instead of the spherical Bessel function. 
Because of this, a correction factor to the diagonal 
elements of 0.5*(diagonal element - 1) is necessary. 

The following text provides the MATLAB program codes 
that I used to calculate the T-matrix. Also provided are 


78 



the user defined functions used in support of the main 
program and the analytical calculation for the T-matrix for 
thin spherical shells. 


Main MATLAB Program : 

% Joseph L. Day 
% Calculation of the T-matrix 
% last updated July 21, 1999 

% Part 1: Import ATILA data and calculate the A coefficients 
format long 

data=getdata('temp5.dat'); 

% terr5>5.dat is a text file, tab deliminated with the first 
% line as; n\imber of rows, one-space, and niomber of columns. 
% The data as follows: colxmin 1: x; col. 2: y; col. 3; z; 

% col. 4: real(pressure); and 
% column 5: Imaginary(pressure). 
x=data(:,1);y=data(:,2);z=data(:,3); 
p=data ( : , 4 ) +i . *data ( : , 5 ) ; 

k=2*pi*474/1492.94; % since ka=l and a=.5meters 

count=l; % initializes my count, 
for N=0:2 

for M=-N:N 

Ptemp2=zeros(1,length(x)); 
for n=l:length(x) 

, r=sqrt((x(n))-2+(y(n))"2+(z(n))''2); 
theta=acos(z(n)/r); 
phi=atan2 (y (n) ,x(n) ) ; 
if M<0 

P= (-1)''abs (M) * fact (N-abs (M)) / 

fact(N+abs(M))*legendre(N,cos(theta)); 

else 

P=legendre(N,cos(theta)); 

end 

Ptempl=shank2(N,k*r)*P(l+abs(M))*exp(i*M*phi); 
Ptemp2(n)=Ptempl; 

end 

Ptemp2=(conj(Ptemp2))'; % make column vector 
eval ([' Ptemp3 ' num2str (coxmt) '=Ptemp2;']); 
count=count+l; 

end 

end 


79 



cal=[Ptemp31 Ptemp32 Ptemp33 Ptemp34 Ptenp35 Ptemp36 Ptemp37 
Ptemp38 Ptemp39]; 

% cal is a 100by9 matrix that will be used to determine 

% the A coefficients 

A=cal\p; 

% This final "A" value provides one column of the T-matrix. 

% At this point I would apply the correction factor since 
% we used a traveling wave instead of the standing wave. 

% The n and m value determines which column of the T-matrix. 


Supporting user defined functions ; 
function A=getdata(name) 

%Function to read free-formated data A (matrix or vector) 
% Written by Charles W. Therrien 
if all(name -='.') 

name=[name,'.dat']; 
end 

[ft,message]=fopen(name,'r'); 
if ft<0 

error(message) 

end 

[bf,N]=fscanf(ft,'%g’); 
fclose(ft); 
if bf(1)==N-1 

A=zeros(l,bf(1)); 

A(:)=bf (2:N) ; 
elseif bf(l)*bf(2)==N-2 
A=zeros(bf(2),bf(1)); 

A(:)=bf(3:N); ' 

A=A. •; 
else 

A=bf.’; 

fprintf([' Data length does not match count.\n',... 

' File contents is being returned as a vector.An']) 

end 


function hn=shank2(n,x) 

% Computes the Spherical Hankel function of the second kind. 

% Written by Joseph L. Day 

hn=sqrt(pi/(2*x))*besselh(n+.5,2,x); 
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function y=fact(n) 

% Computes the factorial of n 
% Written by Joseph L. Day 
y=prod(l:n); 


Program to calculate the Analytical T-matrix for a thin 
spherical shell ; 

% Written by: Joseph L. Day 

% Reference: ONR Report - AY98 Thesis (Ref 5) 

% last updated April 6, 1999 

% Computation of the T-matrix diagonal elements using 
% the analytically derived formulas found (Ref 5, page 5). 
format long 

% is the radius of the sphere in meters 
% is the shell thickness 
% is the bulk modulus in Pascals 
% is the Poisson's ratio shell's material 
% is the angular frequency (474 Hz) 

% is the density of the solid [kg/m'^3] 

% is the density of the fluid [kg/m''3] 

% is the sound speed of the fluid [m/s] 

% is the wave number, in this case ka=l 
cp=sqrt(E/(rhos*(l-v^2))); % plate wave speed 
omega=w*a/cp; 

M=l; 

for N=0:2 
n=N; 

for m=-n:n 

b=omega''2-v-n*(n+l)+l; 
c=omega''2-2* (1+v); 
d=n*(n+l)*(l+v)^2; 

Z(M)=i*h*cp*rhos/(a*omega)*((b*c-d)/b); 
jn=sqrt(pi/(2*k*a))*besselj(n+.5,k*a); 
hn=sqrt(pi/(2*k*a))*besselh(n+.5,2,k*a); 
jhat=sqrt(pi/(2*k*a))*(besselj(n-.5,k*a) 
-(n+l)*besselj(n+.5,k*a)/(k*a)); 
hhat=sqrt(pi/(2*k*a))*(besselh(n-.5,2,k*a) 

-(n+1)*besselh(n+.5,2,k*a)/(k*a)); , 
numerator=i*Z(M)*jhat+rhof*cf*jn; 


a= 0.5; 
h=0.01; 
E=215*10^9; 
■v=0.33; 
w=2*pi*474; 
rhos=7500; 
rhof=1000; 
cf=1492.94; 
k=w/cf; 
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denom=i*Z(M)*hhat+rhof*cf*hn; 

T (M) =-nvraierator/denom; 

M=M+1; 

end 

end 

result=[(1:9)' real(T)' imag(T)' abs(T)* (angle(T)*180/pi)'] 
% result is the output: coliimn 1: index number (n) ; 

% column 2: real(T-matrix); 

% column 3: image(T-matrix); column 4: amplitude; 

% and column 5: phase (degrees). 
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